
> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis ESS
> #Figure 1
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> 
> ####################
> # Load data
> ####################
> source("dataframes.R")

> ##########################
> #Redistribution preferences
> ##########################
> 
> #Average preferences for redistribution (above mean earners)
> #Produces information on redistribution preferences (in percentages) in main text.
> length(data_imputed$redistribution[data_imputed$income.dist.th>0])
[1] 367129

> prefs<-100*prop.table(table(data_imputed$redistribution[data_imputed$income.dist.th>0]))

> rownames(prefs)<-c("Disagree strongly","Disagree","Neither",
+                    "Agree","Agree strongly")

> xtable(t(prefs), caption="Preferences for redistribution (in %)")
% latex table generated in R 4.0.4 by xtable 1.8-4 package
% Wed Aug 10 16:06:17 2022
\begin{table}[ht]
\centering
\begin{tabular}{rrrrrr}
  \hline
 & Disagree strongly & Disagree & Neither & Agree & Agree strongly \\ 
  \hline
1 & 3.82 & 17.12 & 17.94 & 41.16 & 19.96 \\ 
   \hline
\end{tabular}
\caption{Preferences for redistribution (in %)} 
\end{table}

> ##########################
> #Produces Figure 1: Support for redistribution across Western European countries (above mean earners, 2002-2014).
> ##########################
> 
> data_imputed%>%
+   group_by(country,redistribution)%>%
+   filter(income.dist.th>0)%>%
+   select(redistribution)%>%
+   summarise (n = n()) %>%
+   mutate(freq = n / sum(n))->data01

> data01%>%
+   group_by(country)%>%
+   mutate(high=(freq[redistribution==5]+freq[redistribution==4]),
+          low=(freq[redistribution==1]+freq[redistribution==2]),
+          neither=freq[redistribution==3])%>%
+   gather(support, value, c(high,low,neither))->data01

> data01$support[data01$support=="high"]<-3

> data01$support[data01$support=="neither"]<-2

> data01$support[data01$support=="low"]<-1

> data01%>%
+   select(country, support, value)%>%
+   arrange(desc(support),desc(value))%>% 
+   distinct(value,support)->data01

> data01$country <- factor(data01$country, 
+                          levels = unique(data01$country))

> data01%>%
+   ggplot(aes(fill=factor(support),y=value,x=country, 
+              label = paste0(round(value,2)*100,"%")), size=4)+
+   geom_bar(stat = "identity")+
+   geom_text(size = 3, position = position_stack(vjust = 0.5))+
+   coord_flip()+
+   ylab("Redistribution preferences\n in Western European countries (2002-2014).")+
+   theme_bw()+
+   theme(panel.border=element_blank(),axis.line=element_line(),
+         legend.position="bottom",legend.text = element_text(size = 14),
+         legend.title = element_text(size=14),axis.text=element_text(size=12),
+         axis.title=element_text(size=14),
+         axis.text.x = element_blank(),
+         axis.ticks.x = element_blank(),
+         axis.title.y = element_blank())+
+   scale_fill_brewer(palette="Greys",
+                     name=element_blank(),
+                     labels=c("Oppose","Neither nor","Support"),
+                     guide = guide_legend(reverse = TRUE))

> #Compare with listwise deletion
> prefs_list<-100*prop.table(table(data_list$redistribution[data_list$income.dist.th>0]))

> rownames(prefs_list)<-c("Disagree strongly","Disagree","Neither",
+                         "Agree","Agree strongly")

> xtable(t(prefs_list), caption="Preferences for redistribution (in %)")
% latex table generated in R 4.0.4 by xtable 1.8-4 package
% Wed Aug 10 16:06:19 2022
\begin{table}[ht]
\centering
\begin{tabular}{rrrrrr}
  \hline
 & Disagree strongly & Disagree & Neither & Agree & Agree strongly \\ 
  \hline
1 & 3.95 & 18.48 & 17.59 & 40.72 & 19.26 \\ 
   \hline
\end{tabular}
\caption{Preferences for redistribution (in %)} 
\end{table}

> ##########################
> #Save figure
> ##########################
> 
> ggsave(file="figure_1.pdf", height = 5.83, width = 8.27, units = "in")

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis Tax and Benefit System
> #Figure 2
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> 
> 
> ####################
> # Load data
> ####################
> 
> load("Benefits.Rda")

> ##########################
> #Descriptives on benefit concentration indicator
> ##########################
> 
> #Benefit concentration indicator mean
> df %>%
+   group_by(cntry) %>%
+   filter(time=="Mean") %>%
+   summarize(mean = mean(allRepl))
# A tibble: 15 × 2
   cntry   mean
   <chr>  <dbl>
 1 AT    0.385 
 2 BE    0.640 
 3 CH    0.250 
 4 DE    0.0861
 5 DK    0.647 
 6 ES    0.609 
 7 FI    0.368 
 8 FR    0.0932
 9 GB    0.777 
10 IE    0.661 
11 IT    0.450 
12 NL    0.374 
13 NO    0.446 
14 PT    0.138 
15 SE    0.622 

> ##########################
> #Figure 2: Benefit concentration indicator across countries.
> ##########################
> 
> df %>%
+   group_by(cntry) %>%
+   filter(time=="Mean") %>%
+   summarize(mean = mean(allRepl)) %>%
+   arrange(+mean) %>%                               
+   mutate(group = factor(cntry, cntry)) %>%   
+   ggplot(aes(x=group, y=mean)) +                 
+   geom_bar(stat="identity",fill="#333333")+
+   annotate("segment", x = 1.5, xend = 1.5, y = .7, yend = .82,arrow=arrow())+
+   annotate("segment", x = 1.5, xend = 1.5, y = .6, yend = .48,arrow=arrow())+
+   annotate("text", x = c(2.5,2.5), y = c(.75,.55), 
+            label = c("flat-rate", "earnings-related"),hjust=0, size=5) +
+   xlab("Country") + ylab("Benefit concentration") +  coord_cartesian(ylim=c(0,.85))+
+   scale_y_continuous(expand = c(0, 0)) + 
+   theme_bw()+
+   theme(panel.border=element_blank(),axis.line=element_line(),
+         legend.position="bottom",legend.text = element_text(size = 14),
+         legend.title = element_text(size=14),axis.text=element_text(size=12),
+         axis.title=element_text(size=14)
+   )

> ##########################
> 
> ##########################
> #Save figure
> ##########################
> 
> ggsave(file="figure_2.pdf", height = 5.83, width = 8.27, units = "in")

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis ESS
> #Figure 3
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> ####################
> # Load data
> ####################
> source("dataframes.R")

> ##########################
> #Produces Figure 3: Average support for redistribution and benefit concentration (above mean earners, 2002-2014).
> ##########################
> 
> 
> data_imputed_sample%>%
+   filter(income.dist.th>0)%>%
+   ggplot(aes(x = benefitsAvg, y = redistribution,label=country)) +
+   stat_summary(fun.data = mean_se,aes(color="rich"),alpha=.5)+
+   geom_smooth(method=lm,se=F,color="grey20",na.rm=T,size=.5)+
+   xlab("Benefit concentration")+ ylab("Support redistribution")+
+   coord_cartesian(ylim=c(2.5, 4.7), xlim=c(0,.95))+
+   stat_summary(aes(label=country), fun="mean", geom="text", 
+                hjust=0.5, vjust=-.8,size=6)+
+   scale_colour_manual(name="Income group", 
+                       values=c("grey20"),
+                       breaks = c("rich"),
+                       labels = c("rich")
+   )+
+   scale_x_continuous(expand = c(0, 0)) + 
+   scale_y_continuous(expand = c(0, 0)) +
+   theme_bw()+
+   theme(panel.border=element_blank(),axis.line=element_line(),
+         legend.position="none",legend.text = element_text(size = 18),
+         legend.title = element_text(size=18),axis.text=element_text(size=18),
+         axis.title=element_text(size=22)
+   )

> #Check numbers
> data_imputed_sample%>%
+   group_by(country)%>%
+   filter(income.dist.th>=0)%>%
+   summarize(redistribution=mean(redistribution))
# A tibble: 15 × 2
   country        redistribution
   <fct>                   <dbl>
 1 Austria                  3.72
 2 Belgium                  3.55
 3 Denmark                  2.92
 4 Finland                  3.75
 5 France                   3.82
 6 Germany                  3.43
 7 Ireland                  3.80
 8 Italy                    4.02
 9 Netherlands              3.19
10 Norway                   3.37
11 Portugal                 4.20
12 Spain                    3.98
13 Sweden                   3.57
14 Switzerland              3.38
15 United Kingdom           3.30

> ##########################
> #Save figure
> ##########################
> 
> ggsave(file="figure_3.pdf", height = 5.83, width = 8.27, units = "in")

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis ESS
> #Figure A 1
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> 
> ####################
> # Load data
> ####################
> source("dataframes.R")

> ##########################
> #Produces Figure A.1: Distribution of income distance.
> ##########################
> 
> ggplot(data_list, aes(x=income.dist.th)) + 
+   geom_density(kernel="gaussian",bw=14,color="blue") + 
+   geom_density(data = data_imputed, aes(x=income.dist.th), 
+                colour = "red",kernel="gaussian",bw=14) +
+   xlab("Income distance (in 1,000s PPP USD)")+ 
+   ylab("Density")+
+   facet_wrap(~country)+
+   coord_cartesian(xlim=c(-75,125),ylim = c(0,.0225))+
+   theme_bw()+
+   theme(panel.border=element_blank(),axis.line=element_line(),
+         legend.position="bottom",legend.text = element_text(size = 14),
+         legend.title = element_text(size=14),axis.text=element_text(size=12),
+         axis.title=element_text(size=14),strip.background =element_rect(fill=F,colour = "white"))

> ##########################
> #Save figure
> ##########################
> 
> ggsave(file="figure_a1.pdf", height = 5.83, width = 8.27, units = "in")

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis ESS
> #Figure A 2
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> 
> ####################
> # Load data
> ####################
> source("dataframes.R")

> ##########################
> #Produces Figure A.2: Distribution of unemployment risk over 2-digit ISCO.
> ##########################
> 
> 
> summary(data_imputed$isco2d)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -64.42   32.00   51.00   49.95   71.87  157.13 

> ggplot(data_list,aes(x=isco2d, y=our))+
+   xlim(c(0,50))+
+   stat_summary_bin(aes(y = our*100), fun = "mean", 
+                    geom = "bar",fill="blue")+
+   stat_summary_bin(data=data_imputed,aes(y = our*100), 
+                    fun = "mean", geom = "bar",
+                    fill="red",alpha=.5)+
+   xlab("ISCO 2-digit") +  
+   ylab ("Unemployment rate (in %)")+ 
+   theme_bw()+
+   theme(panel.border=element_blank(),axis.line=element_line(),
+         legend.position="bottom",legend.text = element_text(size = 14),
+         legend.title = element_text(size=14),axis.text=element_text(size=12),
+         axis.title=element_text(size=14),strip.background =element_rect(fill=F,colour = "white"))+
+   facet_wrap(~factor(country))

> ##########################
> #Save figure
> ##########################
> 
> ggsave(file="figure_a2.pdf", height = 5.83, width = 8.27, units = "in")

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis Tax and Benefit System
> # Figure B 3
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> 
> 
> ####################
> # Load data
> ####################
> 
> load("DataFile_07_Replacement.Rda")

> load("replacement.Rda")

> ##########################
> #Illustrations: Married couple
> ##########################
> replacement%>%
+   filter(famtype=="married") %>%
+   group_by(Country) %>%
+   summarise(RR_100 = mean(RR[AW==100],na.rm=T),
+             RR_50 = mean(RR[AW==50],na.rm=T),
+             RR_200 = mean(RR[AW==200],na.rm=T)) -> df

> #Rearrange order in ggplot y-axis
> (df %>%
+   group_by(Country) %>%    
+   arrange(+RR_100))$Country -> order 

> df<-reshape(as.data.frame(df), varying=c("RR_50","RR_100","RR_200"), idvar = "group", 
+             direction ="long", sep = "_")

> ##########################
> #Figure B.3: Replacement rates. Married, two children, spouse works 67%, across Western European countries, calculated for 50% AW, AW, 200% AW respectively.
> ##########################
> ggplot(df,aes(RR,Country,color = as.factor(time)))+ 
+   geom_point(size=4,aes(color=as.factor(time),shape=as.factor(time)))+
+   scale_shape_manual(values=c(16,8,16),
+                      name="AW (in %)",
+                      guide = guide_legend(reverse = TRUE))+
+   scale_color_manual(values=c("lightgray","black","darkgray"),
+                      name="AW (in %)",
+                      guide = guide_legend(reverse = TRUE))+
+   scale_y_discrete(limits = order)+
+   ylab("Country") + xlab("Replacement rate") +
+   theme_bw()+
+   theme(panel.border=element_blank(),
+         axis.line=element_line(),
+         legend.position="bottom",
+         legend.text = element_text(size = 14),
+         legend.title = element_text(size=14),
+         axis.text=element_text(size=12),
+         axis.title=element_text(size=14))

> ##########################
> #Save figure
> ##########################
> 
> ggsave(file="figure_b3.pdf", height = 5.83, width = 8.27, units = "in")

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis Tax and Benefit System
> # Figure B 4
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> 
> 
> ####################
> # Load data
> ####################
> 
> load("DataFile_07_Replacement.Rda")

> load("replacement.Rda")

> ##########################
> #Illustrations: Single earner
> ##########################
> 
> replacement%>%
+   filter(famtype=="single") %>%
+   group_by(Country) %>%
+   summarise(RR_100 = mean(RR[AW==100],na.rm=T),
+             RR_50 = mean(RR[AW==50],na.rm=T),
+             RR_200 = mean(RR[AW==200],na.rm=T)) -> df

> #Rearrange order in ggplot y-axis
> (df %>%
+     group_by(Country) %>%    
+     arrange(+RR_100))$Country -> order 

> df<-reshape(as.data.frame(df), varying=c("RR_50","RR_100","RR_200"), idvar = "group", 
+             direction ="long", sep = "_")

> ##########################
> #Figure B.4: Replacement rates. Single, no children, across Western European countries, calculated for 50% AW, AW, 200% AW respectively.
> ##########################
> ggplot(df,aes(RR,Country,color = as.factor(time)))+ 
+   geom_point(size=4,aes(color=as.factor(time),shape=as.factor(time)))+
+   scale_shape_manual(values=c(16,8,16),
+                      name="AW (in %)",
+                      guide = guide_legend(reverse = TRUE))+
+   scale_color_manual(values=c("lightgray","black","darkgray"),
+                      name="AW (in %)",
+                      guide = guide_legend(reverse = TRUE))+
+   scale_y_discrete(limits = order)+
+   ylab("Country") + xlab("Replacement rate") +
+   theme_bw()+
+   theme(panel.border=element_blank(),
+         axis.line=element_line(),
+         legend.position="bottom",
+         legend.text = element_text(size = 14),
+         legend.title = element_text(size=14),
+         axis.text=element_text(size=12),
+         axis.title=element_text(size=14))

> ##########################
> #Save figure
> ##########################
> 
> ggsave(file="figure_b4.pdf", height = 5.83, width = 8.27, units = "in")

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis Tax and Benefit System
> # Figure B 5
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> 
> 
> ####################
> # Load data
> ####################
> load("Benefits.Rda")

> #Combine benefit concentration indicator for each family type in one plot
> #Order by mean replacement rate
> (df %>%
+    group_by(cntry) %>% 
+    filter(time=="Mean") %>%
+    summarize(mean = mean(allRepl)) %>%
+    arrange(+mean))$cntry -> order

> #Produces A.2 in Appendix
> ggplot(df,aes(allRepl,cntry,color = as.factor(time)))+ 
+   geom_point(size=4,aes(shape=as.factor(time)))+
+   scale_color_manual(values=c("lightgray","black","darkgray"),
+                      name="Family type")+
+   scale_shape_manual(values=c(16,8,16),name="Family type")+
+   scale_y_discrete(limits = order)+
+   ylab("Country") + xlab("Benefit concentration indicator") +
+   theme_bw()+
+   theme(panel.border=element_blank(),
+         axis.line=element_line(),
+         legend.position="bottom",
+         legend.text = element_text(size = 14),
+         legend.title = element_text(size=14),
+         axis.text=element_text(size=12),
+         axis.title=element_text(size=14))

> ##########################
> #Save figure
> ##########################
> 
> ggsave(file="figure_b5.pdf", height = 5.83, width = 8.27, units = "in")

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis ESS
> #Figure B 6
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> ####################
> # Load data
> ####################
> source("dataframes.R")

> ##########################
> #Benefit concentration indicators versus alternative measures: Benefit concentration and welfare spending
> ##########################
> 
> cor.test(data_sub$social_spending,data_sub$benefitsAvgYear, method=c("pearson"))

	Pearson's product-moment correlation

data:  data_sub$social_spending and data_sub$benefitsAvgYear
t = -88.182, df = 192366, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2014028 -0.1928127
sample estimates:
       cor 
-0.1971115 


> ##########################
> #Figure B.6: Benefit concentration and social spending (% of GDP), averaged over 2002- 2014.
> ##########################
> 
> data_sub%>%
+   group_by(country)%>%
+   summarize(benefits=mean(benefitsAvg),
+             social_spending=mean(social_spending))%>%
+   ggplot(aes(x=benefits, y=social_spending)) + 
+   geom_point() + 
+   #ggtitle("Benefit concentration and social expenditure (total, % of GDP)") + 
+   geom_smooth(method=lm, se=T) + 
+   scale_x_continuous(name = "Benefit concentration", limits = c(0, 1)) + 
+   scale_y_continuous(name = "Social spending (% of GDP)") + 
+   theme_bw()

> ##########################
> #Save figure
> ##########################
> 
> ggsave(file="figure_b6.pdf", height = 5.83, width = 8.27, units = "in")

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis ESS
> #Figure B 7
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> ####################
> # Load data
> ####################
> source("dataframes.R")

> ##########################
> #Benefit concentration indicators versus alternative measures: Benefit concentration and top statutory personal income tax rate
> ##########################
> 
> cor.test(data_sub$social_spending,data_sub$benefitsAvgYear, method=c("pearson"))

	Pearson's product-moment correlation

data:  data_sub$social_spending and data_sub$benefitsAvgYear
t = -88.182, df = 192366, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2014028 -0.1928127
sample estimates:
       cor 
-0.1971115 


> ##########################
> # Produces Figure B.7: Benefit concentration and top statutory personal income tax rate, averaged over 2002-2014.
> ##########################
> 
> data_sub%>%
+   group_by(country)%>%
+   summarize(benefits=mean(benefitsAvg),
+             toptax=mean(toptax))%>%
+   ggplot(aes(x=benefits, y=toptax)) +
+   geom_point() + 
+   #ggtitle("Benefit concentration and top income tax") + 
+   geom_smooth(method=lm, se=T) + 
+   scale_x_continuous(name = "Benefit concentration", limits = c(0, 1)) + 
+   scale_y_continuous(name = "Top income tax") + 
+   theme_bw()

> ##########################
> #Save figure
> ##########################
> 
> ggsave(file="figure_b7.pdf", height = 5.83, width = 8.27, units = "in")

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis ESS
> #Figure B 8
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> ####################
> # Load data
> ####################
> source("dataframes.R")

> ##########################
> #Benefit concentration indicators versus alternative measures: Benefit concentration and fiscal redistribution (benefits) by Mahler and Jesuit (2006)
> ##########################
> 
> cor.test(data_sub$red_transfers,data_sub$benefitsAvgYear, method=c("pearson"))

	Pearson's product-moment correlation

data:  data_sub$red_transfers and data_sub$benefitsAvgYear
t = 38.844, df = 38912, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1836254 0.2027551
sample estimates:
      cor 
0.1932086 


> ##########################
> # Produces Figure B.8: Benefit concentration and fiscal redistribution achieved via transfers.
> ##########################
> 
> 
> data_sub%>%
+   group_by(country)%>%
+   summarize(benefits=mean(benefitsAvg),
+             red_transfers=mean(red_transfers,na.rm=T))%>%
+   ggplot(aes(x=benefits, y=red_transfers)) + 
+   geom_point() + 
+   #ggtitle("Benefit concentration and fiscal redistribution (via transfers)") + 
+   geom_smooth(method=lm, se=T) + 
+   scale_x_continuous(name = "Benefit concentration", limits = c(0, 1)) + 
+   scale_y_continuous(name = "Fiscal redistribution (via transfers)") + 
+   theme_bw()

> ##########################
> #Save figure
> ##########################
> 
> ggsave(file="figure_b8.pdf", height = 5.83, width = 8.27, units = "in")

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis ESS
> #Figure B 9
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> ####################
> # Load data
> ####################
> source("dataframes.R")

> ##########################
> #Benefit concentration indicators versus alternative measures: Benefit concentration and target efficiency by Mahler and Jesuit (2006)
> ##########################
> data_sub$efficiency<-NA

> data_sub$efficiency[data_sub$country=="Belgium"]<- -0.125

> data_sub$efficiency[data_sub$country=="Sweden"]<- -0.062

> data_sub$efficiency[data_sub$country=="Netherlands"]<- -0.015

> data_sub$efficiency[data_sub$country=="Finland"]<- -0.134

> data_sub$efficiency[data_sub$country=="France"]<- 0.053

> data_sub$efficiency[data_sub$country=="Denmark"]<- -0.142

> data_sub$efficiency[data_sub$country=="Germany"]<- -0.223

> data_sub$efficiency[data_sub$country=="United Kingdom"]<- -0.281

> data_sub$efficiency[data_sub$country=="Norway"]<- -0.246

> data_sub$efficiency[data_sub$country=="Switzerland"]<- -0.059

> cor.test(data_sub$efficiency,data_sub$benefitsAvg, method=c("pearson"))

	Pearson's product-moment correlation

data:  data_sub$efficiency and data_sub$benefitsAvg
t = -126.28, df = 137206, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.3274142 -0.3179335
sample estimates:
       cor 
-0.3226819 


> ##########################
> # Produces Figure B.9: Benefit concentration and internal target efficiency of social benefits.
> ##########################
> 
> 
> data_sub%>%
+   group_by(country)%>%
+   filter(year==2002)%>%
+   summarize(benefits=mean(benefitsAvg),
+             efficiency=mean(efficiency,na.rm=T))%>%
+   ggplot(aes(x=benefits, y=efficiency)) + 
+   geom_point() + 
+   ggtitle("Benefit concentration and fiscal redistribution (via transfers)") + 
+   geom_smooth(method=lm, se=T) + 
+   scale_x_continuous(name = "Benefit concentration", limits = c(0, 1)) + 
+   scale_y_continuous(name = "Target efficiency") + 
+   theme_bw()

> ##########################
> #Save figure
> ##########################
> 
> ggsave(file="figure_b9.pdf", height = 5.83, width = 8.27, units = "in")

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis Tax and Benefit System
> #Table B 2
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> 
> 
> ####################
> # Load data
> ####################
> 
> load("DataFile_07_Replacement.Rda")

> load("replacement.Rda")

> ##########################
> #Illustrations: Married couple
> ##########################
> 
> replacement%>%
+   filter(famtype=="married") %>%
+   group_by(Country) %>%
+   summarise(RR_100 = mean(RR[AW==100],na.rm=T),
+             RR_50 = mean(RR[AW==50],na.rm=T),
+             RR_200 = mean(RR[AW==200],na.rm=T)) -> df

> min<-rbind(min(df$RR_50),min(df$RR_100),min(df$RR_200))

> mean<-rbind(mean(df$RR_50),mean(df$RR_100),mean(df$RR_200))

> max<-rbind(max(df$RR_50),max(df$RR_100),max(df$RR_200))

> df_table<-cbind(min,mean,max)

> rownames(df_table)<-c("AW 50%", "AW 100%", "AW 200%")

> colnames(df_table)<-c("min", "mean", "max")

> ##########################
> #Produces Table B.2: Married couple, 2 children, gross earning spouse 67% of AW.
> ##########################
> xtable(df_table, caption="Married couple, 2 children, gross earning spouse 67% of AW")
% latex table generated in R 4.0.4 by xtable 1.8-4 package
% Wed Aug 10 16:06:51 2022
\begin{table}[ht]
\centering
\begin{tabular}{rrrr}
  \hline
 & min & mean & max \\ 
  \hline
AW 50\% & 0.75 & 0.88 & 0.98 \\ 
  AW 100\% & 0.56 & 0.79 & 0.92 \\ 
  AW 200\% & 0.39 & 0.60 & 0.80 \\ 
   \hline
\end{tabular}
\caption{Married couple, 2 children, gross earning spouse 67% of AW} 
\end{table}

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis: ESS
> # Table B 3
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> 
> 
> ####################
> # Load data
> ####################
> 
> load("DataFile_07_Replacement.Rda")

> load("replacement.Rda")

> ##########################
> #Illustrations: Single earner
> ##########################
> 
> 
> replacement%>%
+   filter(famtype=="single") %>%
+   group_by(Country) %>%
+   summarise(RR_100 = mean(RR[AW==100],na.rm=T),
+             RR_50 = mean(RR[AW==50],na.rm=T),
+             RR_200 = mean(RR[AW==200],na.rm=T)) -> df

> min<-rbind(min(df$RR_50),min(df$RR_100),min(df$RR_200))

> mean<-rbind(mean(df$RR_50),mean(df$RR_100),mean(df$RR_200))

> max<-rbind(max(df$RR_50),max(df$RR_100),max(df$RR_200))

> df_table<-cbind(min,mean,max)

> rownames(df_table)<-c("AW 50%", "AW 100%", "AW 200%")

> colnames(df_table)<-c("min", "mean", "max")

> ##########################
> #Produces Table B.3: Single person, no children.
> ##########################
> xtable(df_table, caption="Single person, no children.")
% latex table generated in R 4.0.4 by xtable 1.8-4 package
% Wed Aug 10 16:06:51 2022
\begin{table}[ht]
\centering
\begin{tabular}{rrrr}
  \hline
 & min & mean & max \\ 
  \hline
AW 50\% & 0.66 & 0.78 & 0.94 \\ 
  AW 100\% & 0.39 & 0.61 & 0.80 \\ 
  AW 200\% & 0.22 & 0.42 & 0.68 \\ 
   \hline
\end{tabular}
\caption{Single person, no children.} 
\end{table}

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis ESS
> #PTable C 4
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> 
> ####################
> # Load data
> ####################
> source("dataframes.R")

> ##########################
> #Summary statistics
> ##########################
> summean<-NULL

> for(i in 1:5) {
+   sumry <- mi.out[[i]] %>%
+     summarise(red=mean(redistribution),
+               inc=mean(income.PPP.th),
+               incdist=mean(income.dist.th),
+               our=mean(our),
+               age=mean(agea),
+               education=mean(eduyrs),
+               benefit=mean(benefitsAvg),
+               gini=mean(gini.disp),
+               immigration=mean(immigration),
+               social_spending=mean(social_spending),
+               toptax=mean(toptax),
+               red_transfers=mean(red_transfers,na.rm=T),
+               single=mean(single)*100,
+               gender=mean(gender)*100,
+               unemployment=mean(uempl)*100)
+   summean <- rbind(summean, sumry[,1:15])
+ }

> sumsd<-NULL

> for(i in 1:5) {
+   sumry <- mi.out[[i]] %>%
+     summarise(redsd=sd(redistribution),
+               incsd=sd(income.PPP.th),
+               incdistsd=sd(income.dist.th),
+               oursd=sd(our),
+               agesd=sd(agea),
+               educationsd=sd(eduyrs),
+               benefitsd=sd(benefitsAvg),
+               ginisd=sd(gini.disp),
+               immigrationsd=sd(immigration),
+               social_spending=sd(social_spending),
+               toptax=sd(toptax),
+               red_transfers=sd(red_transfers,na.rm=T),
+               singlesd=0,
+               gendersd=0,
+               unemploymentsd=0)
+   sumsd <- rbind(sumsd, sumry[,1:15])
+ }

> summean[3]
        incdist
1  1.228143e-17
2 -1.001396e-16
3  2.067626e-16
4 -1.672403e-16
5  1.522508e-16

> sumsd[3]
  incdistsd
1  27.18376
2  27.22816
3  27.17571
4  27.13502
5  27.17913

> allmean<-NULL

> allsd<-NULL

> for(i in 1:15) {
+   combined.results <- mi.meld(q = summean[i], se = sumsd[i])
+   allmean <- rbind(allmean,combined.results$q.mi)
+   allsd <- rbind(allsd,combined.results$se.mi)
+ }

> df_imp<-data.frame(cbind(round(allmean,2),round(allsd,2)))

> data_list %>%
+   summarise(red=mean(redistribution),
+             inc=mean(income.PPP.th),
+             incdist=mean(income.dist.th),
+             our=mean(our),
+             age=mean(agea),
+             education=mean(eduyrs),
+             benefit=mean(benefitsAvg),
+             gini=mean(gini.disp),
+             immigration=mean(immigration),
+             social_spending=mean(social_spending),
+             toptax=mean(toptax),
+             red_transfers=mean(red_transfers,na.rm=T),
+             single=mean(single)*100,
+             gender=mean(gender)*100,
+             unemployment=mean(uempl)*100) -> summean_list

> data_list %>%
+   summarise(red=sd(redistribution),
+             inc=sd(income.PPP.th),
+             incdist=sd(income.dist.th),
+             our=sd(our),
+             age=sd(agea),
+             education=sd(eduyrs),
+             benefit=sd(benefitsAvg),
+             gini=sd(gini.disp),
+             immigration=sd(immigration),
+             social_spending=sd(social_spending),
+             toptax=sd(toptax),
+             red_transfers=sd(red_transfers,na.rm=T),
+             single=0,
+             gender=0,
+             unemployment=0) -> sumsd_list

> df_list<-t(data.frame(rbind(round(summean_list,2),round(sumsd_list,2))))

> df<-cbind(df_imp,df_list)

> rownames(df)<-c("Redistribution","Income [1,000 PPP USD]",
+                 "Income distance","Occupational unemployment [in %]",
+                 "Age [in years]", "Education [in years]",
+                 "Benefit concentration","Gini(disposable income)",
+                 "Immigration","Social spending", "Top tax","Redistribution via transfers",
+                 "Single","Female","Unemployed")

> colnames(df)<-c("Mean","SD","Mean","SD")

> ##########################
> #Produces Table C.4: Descriptive statistics. Means, standard deviations, and percentages.
> ##########################
> 
> xtable(df, caption="Descriptive statistics (comparing multiple 
+ imputation and listwise deletion). Means, standard deviations 
+        and percentages.")
% latex table generated in R 4.0.4 by xtable 1.8-4 package
% Wed Aug 10 16:06:54 2022
\begin{table}[ht]
\centering
\begin{tabular}{rrrrr}
  \hline
 & Mean & SD & Mean & SD \\ 
  \hline
Redistribution & 3.75 & 1.05 & 3.73 & 1.06 \\ 
  Income [1,000 PPP USD] & 36.69 & 28.34 & 37.41 & 28.26 \\ 
  Income distance & 0.00 & 27.18 & 0.52 & 26.97 \\ 
  Occupational unemployment [in \%] & 0.07 & 0.06 & 0.06 & 0.06 \\ 
  Age [in years] & 48.16 & 18.52 & 49.01 & 17.11 \\ 
  Education [in years] & 12.41 & 4.30 & 12.75 & 4.20 \\ 
  Benefit concentration & 0.43 & 0.23 & 0.43 & 0.23 \\ 
  Gini(disposable income) & 28.59 & 3.42 & 28.17 & 3.30 \\ 
  Immigration & 11.69 & 4.93 & 11.56 & 4.94 \\ 
  Social spending & 49.99 & 4.39 & 49.82 & 4.42 \\ 
  Top tax & 48.50 & 5.53 & 48.75 & 5.67 \\ 
  Redistribution via transfers & 0.16 & 0.02 & 0.17 & 0.02 \\ 
  Single & 22.54 & 0.00 & 18.75 & 0.00 \\ 
  Female & 52.36 & 0.00 & 50.63 & 0.00 \\ 
  Unemployed & 4.16 & 0.00 & 4.03 & 0.00 \\ 
   \hline
\end{tabular}
\caption{Descriptive statistics (comparing multiple 
imputation and listwise deletion). Means, standard deviations 
       and percentages.} 
\end{table}

> ################
> #PSRM: Explaining Support for Redistribution: Social Insurance Systems and Fairness
> #
> #Observational Analysis ESS
> #Table C 5
> #
> #Verena Fetscher
> #July 2022
> ####################
> 
> ####################
> # Load data
> ####################
> source("dataframes.R")

> ##########################
> #Model estimation
> ##########################
> 
> ##########################
> #Fixed effects
> ##########################
> mod1<-zelig(redistribution ~ benefitsAvg+
+               factor(year),
+             model = 'ls',
+             data=mi.out$imp5[mi.out$imp5$income.dist.th>0,], 
+             cite = FALSE)

> summary(mod1)
Model: 

Call:
z5$zelig(formula = redistribution ~ benefitsAvg + factor(year), 
    data = mi.out$imp5[mi.out$imp5$income.dist.th > 0, ])

Residuals:
    Min      1Q  Median      3Q     Max 
-2.7606 -0.6401  0.3615  0.5410  1.6322 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.62512    0.01372 264.216  < 2e-16
benefitsAvg -0.33130    0.01789 -18.519  < 2e-16
year2004     0.04602    0.01537   2.995  0.00275
year2006     0.06275    0.01542   4.070 4.71e-05
year2008     0.08026    0.01554   5.165 2.41e-07
year2010     0.07059    0.01596   4.424 9.72e-06
year2012     0.16403    0.01531  10.710  < 2e-16
year2014     0.13719    0.01549   8.855  < 2e-16

Residual standard error: 1.099 on 73357 degrees of freedom
Multiple R-squared:  0.006783,	Adjusted R-squared:  0.006688 
F-statistic: 71.57 on 7 and 73357 DF,  p-value: < 2.2e-16

Next step: Use 'setx' method

> mod2<-zelig(redistribution ~ benefitsAvg+
+               our+agea+eduyrs+single+gender+uempl+
+               factor(year),
+             model = 'ls',
+             data=mi.out$imp5[mi.out$imp5$income.dist.th>0,],
+             cite = FALSE)

> summary(mod2)
Model: 

Call:
z5$zelig(formula = redistribution ~ benefitsAvg + our + agea + 
    eduyrs + single + gender + uempl + factor(year), data = mi.out$imp5[mi.out$imp5$income.dist.th > 
    0, ])

Residuals:
    Min      1Q  Median      3Q     Max 
-4.8822 -0.6842  0.2907  0.6814  2.0361 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.3387986  0.0267108 124.998  < 2e-16
benefitsAvg -0.2446534  0.0176247 -13.881  < 2e-16
our          3.4019396  0.0833955  40.793  < 2e-16
agea         0.0030638  0.0002991  10.244  < 2e-16
eduyrs      -0.0149740  0.0010156 -14.744  < 2e-16
single       0.0513007  0.0117891   4.352 1.35e-05
gender       0.2059862  0.0079770  25.822  < 2e-16
uempl        0.1762174  0.0280652   6.279 3.43e-10
year2004     0.0346810  0.0150368   2.306   0.0211
year2006     0.0713608  0.0150965   4.727 2.28e-06
year2008     0.0847958  0.0152186   5.572 2.53e-08
year2010     0.0869526  0.0156523   5.555 2.78e-08
year2012     0.1606599  0.0150505  10.675  < 2e-16
year2014     0.1409000  0.0152503   9.239  < 2e-16

Residual standard error: 1.075 on 73351 degrees of freedom
Multiple R-squared:  0.04945,	Adjusted R-squared:  0.04929 
F-statistic: 293.6 on 13 and 73351 DF,  p-value: < 2.2e-16

Next step: Use 'setx' method

> mod3<-zelig(redistribution ~ benefitsAvg+
+               our+agea+eduyrs+single+gender+uempl+
+               gini.disp+immigration+social_spending+toptax+red_transfers+
+               as.factor(year),
+             model = 'ls',
+             data=mi.out$imp5[mi.out$imp5$income.dist.th>0,],
+             cite = FALSE)

> summary(mod3)
Model: 

Call:
z5$zelig(formula = redistribution ~ benefitsAvg + our + agea + 
    eduyrs + single + gender + uempl + gini.disp + immigration + 
    social_spending + toptax + red_transfers + as.factor(year), 
    data = mi.out$imp5[mi.out$imp5$income.dist.th > 0, ])

Residuals:
    Min      1Q  Median      3Q     Max 
-3.0895 -0.7982  0.2435  0.7269  2.4575 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)
(Intercept)      4.6285371  0.2286085  20.247  < 2e-16
benefitsAvg     -0.4240487  0.0566227  -7.489 7.32e-14
our              2.9901863  0.2132852  14.020  < 2e-16
agea             0.0022589  0.0006797   3.323 0.000892
eduyrs          -0.0176447  0.0023885  -7.387 1.58e-13
single           0.0798088  0.0259520   3.075 0.002107
gender           0.2132881  0.0175993  12.119  < 2e-16
uempl            0.1991303  0.0598914   3.325 0.000887
gini.disp        0.0396959  0.0040803   9.729  < 2e-16
immigration     -0.0179637  0.0022924  -7.836 4.95e-15
social_spending -0.0298150  0.0021346 -13.967  < 2e-16
toptax          -0.0325191  0.0025436 -12.785  < 2e-16
red_transfers    5.6794552  0.5613549  10.117  < 2e-16
year2004         0.0008442  0.0512680   0.016 0.986862
year2010        -0.0067405  0.0559583  -0.120 0.904123

Residual standard error: 1.073 on 15118 degrees of freedom
  (58232 observations deleted due to missingness)
Multiple R-squared:  0.07255,	Adjusted R-squared:  0.07169 
F-statistic: 84.47 on 14 and 15118 DF,  p-value: < 2.2e-16

Next step: Use 'setx' method

> texreg(list(mod1,mod2,mod3),stars = c(0.01,0.05,0.1),
+        custom.coef.names = c("Intercept", "Benefit concentration","Year 2004",
+                              "Year 2006","Year 2008","Year 2010","Year 2012",
+                              "Year 2014","Occu. unemployment risk","Age","Education","Single","Female","Unemployment",
+                              "Gini (disposable income)", "Immigration",
+                              "Social Spending","Top Income Tax",
+                              "Redistribution via transfers"))

\begin{table}
\begin{center}
\begin{tabular}{l c c c}
\hline
 & Model 1 & Model 2 & Model 3 \\
\hline
Intercept                    & $3.63^{***}$  & $3.34^{***}$  & $4.63^{***}$  \\
                             & $(0.01)$      & $(0.03)$      & $(0.23)$      \\
Benefit concentration        & $-0.33^{***}$ & $-0.24^{***}$ & $-0.42^{***}$ \\
                             & $(0.02)$      & $(0.02)$      & $(0.06)$      \\
Year 2004                    & $0.05^{***}$  & $0.03^{**}$   & $0.00$        \\
                             & $(0.02)$      & $(0.02)$      & $(0.05)$      \\
Year 2006                    & $0.06^{***}$  & $0.07^{***}$  &               \\
                             & $(0.02)$      & $(0.02)$      &               \\
Year 2008                    & $0.08^{***}$  & $0.08^{***}$  &               \\
                             & $(0.02)$      & $(0.02)$      &               \\
Year 2010                    & $0.07^{***}$  & $0.09^{***}$  & $-0.01$       \\
                             & $(0.02)$      & $(0.02)$      & $(0.06)$      \\
Year 2012                    & $0.16^{***}$  & $0.16^{***}$  &               \\
                             & $(0.02)$      & $(0.02)$      &               \\
Year 2014                    & $0.14^{***}$  & $0.14^{***}$  &               \\
                             & $(0.02)$      & $(0.02)$      &               \\
Occu. unemployment risk      &               & $3.40^{***}$  & $2.99^{***}$  \\
                             &               & $(0.08)$      & $(0.21)$      \\
Age                          &               & $0.00^{***}$  & $0.00^{***}$  \\
                             &               & $(0.00)$      & $(0.00)$      \\
Education                    &               & $-0.01^{***}$ & $-0.02^{***}$ \\
                             &               & $(0.00)$      & $(0.00)$      \\
Single                       &               & $0.05^{***}$  & $0.08^{***}$  \\
                             &               & $(0.01)$      & $(0.03)$      \\
Female                       &               & $0.21^{***}$  & $0.21^{***}$  \\
                             &               & $(0.01)$      & $(0.02)$      \\
Unemployment                 &               & $0.18^{***}$  & $0.20^{***}$  \\
                             &               & $(0.03)$      & $(0.06)$      \\
Gini (disposable income)     &               &               & $0.04^{***}$  \\
                             &               &               & $(0.00)$      \\
Immigration                  &               &               & $-0.02^{***}$ \\
                             &               &               & $(0.00)$      \\
Social Spending              &               &               & $-0.03^{***}$ \\
                             &               &               & $(0.00)$      \\
Top Income Tax               &               &               & $-0.03^{***}$ \\
                             &               &               & $(0.00)$      \\
Redistribution via transfers &               &               & $5.68^{***}$  \\
                             &               &               & $(0.56)$      \\
\hline
R$^2$                        & $0.01$        & $0.05$        & $0.07$        \\
Adj. R$^2$                   & $0.01$        & $0.05$        & $0.07$        \\
Num. obs.                    & $73365$       & $73365$       & $15133$       \\
\hline
\multicolumn{4}{l}{\scriptsize{$^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.1$}}
\end{tabular}
\caption{Statistical models}
\label{table:coefficients}
\end{center}
\end{table}

> ##########################
> #Quantities of interest
> ##########################
> 
> # min to max
> mi.out$imp5%>%
+   group_by(country)%>%
+   summarize(benefits=mean(benefitsAvg))
# A tibble: 15 × 2
   country        benefits
   <fct>             <dbl>
 1 Austria          0.385 
 2 Belgium          0.640 
 3 Denmark          0.647 
 4 Finland          0.368 
 5 France           0.0932
 6 Germany          0.0861
 7 Ireland          0.661 
 8 Italy            0.450 
 9 Netherlands      0.374 
10 Norway           0.446 
11 Portugal         0.138 
12 Spain            0.609 
13 Sweden           0.622 
14 Switzerland      0.250 
15 United Kingdom   0.777 

> x_Germany <- setx(mod2, benefitsAvg=0.0861) 

> x_UK <- setx(mod2, benefitsAvg=0.777) 

> s_out <- sim(mod2, x = x_Germany, x1 = x_UK)

> summary(s_out)

 sim x :
 -----
ev
      mean         sd      50%     2.5%    97.5%
1 3.725005 0.01160687 3.724936 3.701835 3.748786
pv
         mean       sd      50%     2.5%    97.5%
[1,] 3.752556 1.087908 3.725694 1.735766 5.918872

 sim x1 :
 -----
ev
      mean         sd      50%     2.5%    97.5%
1 3.556079 0.01194728 3.555842 3.534554 3.580627
pv
         mean       sd      50%     2.5%    97.5%
[1,] 3.563672 1.109934 3.519124 1.360029 5.765812
fd
        mean         sd        50%       2.5%      97.5%
1 -0.1689256 0.01258685 -0.1692255 -0.1927741 -0.1437625

> # mean to mean+sd
> mi.out$imp5%>%
+   summarize(mean_benefits=mean(benefitsAvg),
+             sd_benefits=sd(benefitsAvg)) -> benefits

> benefits
  mean_benefits sd_benefits
1     0.4283177    0.228922

> benefits[[2]]
[1] 0.228922

> x_mean <- setx(mod2, benefitsAvg=benefits[[1]]) 

> x_sd <- setx(mod2, benefitsAvg=(benefits[[1]]+benefits[[2]]))

> s_out <- sim(mod2, x = x_mean, x1 = x_sd)

> summary(s_out)

 sim x :
 -----
ev
      mean         sd      50%     2.5%    97.5%
1 3.642006 0.01039712 3.641914 3.620977 3.662838
pv
         mean       sd     50%     2.5%    97.5%
[1,] 3.654974 1.090966 3.68776 1.542866 5.757242

 sim x1 :
 -----
ev
      mean         sd      50%     2.5%    97.5%
1 3.585851 0.01127151 3.586067 3.562513 3.606962
pv
         mean       sd      50%     2.5%    97.5%
[1,] 3.568495 1.041847 3.574946 1.506271 5.562447
fd
         mean          sd         50%        2.5%       97.5%
1 -0.05615539 0.004088287 -0.05617133 -0.06470136 -0.04852108
